#include "scrub.hpp"
#include "lencode.hpp"
#include "ldecode.hpp"
#include "bencode.hpp"
#include "bdecode.hpp"
#include "huffman.hpp"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#define _USE_MATH_DEFINES  // in order to use PI
#include <math.h>
#include <set>

namespace Algorithms{

  // for preventing the seed to be the same for each gsl
  time_t gsl_seed = time(NULL) * rand();

  // construct one random number generator
  gsl_rng* get_one_rand_num_generator() {
    const gsl_rng_type * T;
    gsl_rng * r;
    /* create a generator chosen by the 
       environment variable GSL_RNG_TYPE */
    gsl_rng_env_setup();
    gsl_rng_default_seed =  gsl_seed;
    gsl_seed += rand(); // ok for PC, not for HPC: time(NULL) + N * clock() % K;
    sleep(0.2);
#ifdef VERBOSE
    cout << "seed of gsl:" << gsl_rng_default_seed << endl;
#endif
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
    return r;
  }

  // construct a vector of random number generator
  vector<gsl_rng*>* get_rand_num_generators() {
    vector<gsl_rng*>* result = new vector<gsl_rng*>;
    for(int i = 1; i <= MAX_CELL_LEVEL + 1; i++){
      result -> push_back(get_one_rand_num_generator());
    }
    return result;
  }

  // random number generators for R0 and V of each cell level
  // the index starts from 1
  vector<gsl_rng*>* randvarsR = get_rand_num_generators();
  vector<gsl_rng*>* randvarsV = get_rand_num_generators();

  // generate one new initial resistance for the input cell level 
  inline double get_new_resistance(int& level){
    double r = (resistance_means[level] + gsl_ran_gaussian((*randvarsR)[level], resistance_dev));
#ifdef VERBOSE
    if ((level > 1) && ( level < 4))
      if ((r < resistance_thresholds[level - 2]) || (r > resistance_thresholds[level - 1]))
	cout << "Warning: resistance is not precise: " << r << " " << level << endl;
#endif
    return r;
  }

  // generate one new drift exponent for the input cell level 
  inline double get_new_drift_exponent(int& level){
    return (exponent_means[level] + gsl_ran_gaussian((*randvarsV)[level], exponent_devs[level]));
  }

  // construct a new cell with initial resistance and drift exponent
  // corresponding to its level.
  inline Cell* init_one_cell(int& level, vec_GF2* bits){
    double R0 = resistance_means[level] + gsl_ran_gaussian((*randvarsR)[level], resistance_dev);
    double V  = exponent_means[level] + gsl_ran_gaussian((*randvarsV)[level], exponent_devs[level]);
    int L0 = level;
    return (new Cell(R0, L0, V, bits));
  }

  // get the corresponding symbols for the (i+1)-th cell
  inline vec_GF2* get_corresponding_bits(const vec_GF2& codeword, long& i){
    vec_GF2* result = new vec_GF2();
    result -> SetLength(BITS_PER_CELL);
    for(long j = i * BITS_PER_CELL, k = 1; j <= (i + 1) * BITS_PER_CELL - 1; j++, k++){
      (*result)(k) = codeword(j+1);
    } 
    return result;
  }

  // initialize the cell vector to contain the codeword
  void initialize_PCM_cells(vector<Cell*>& cells, const vec_GF2& codeword){
    vector<int> levels;            // index starts from 0
    vector<double> resistances;
    vector<double> exponents;
    // compute cell levels
    CellUtil::compute_cell_levels(levels, codeword);
    for(long i = 0; i < (long) levels.size(); i++){
      vec_GF2* bits = get_corresponding_bits(codeword, i);
      cells.push_back(init_one_cell(levels[i], bits));
      //delete bits;
    }
  }

  void free_PCM_cells(vector<Cell*>& cells){
    for(long i = 0; i < (long) cells.size(); i++){
      delete cells[i] -> bits;
      delete cells[i] -> bits0;
      delete cells[i];
    }
  }

  // compute the value of drifted resistance for a cell at its life time t
  inline double get_drifted_resistance(const Cell* cell){
    return (cell -> R0 + (cell -> V) * log10((double) cell -> t));
  }

  // compute the new level corresponding to the drifted resistance r
  inline int compute_new_level(const double& r){
    int result = 1;
    for(int i = 0; i < (MAX_CELL_LEVEL - 1); i++, result++){
      if (r < resistance_thresholds[i])
	return result;
    }
    return result;
  }

  // compute the stored bits according to the current cell level.
  void update_bits_in_cell(Cell* cell){
    int l = cell -> L;
    if (l == 1) {
      (*(cell -> bits))(1) = 0;
      (*(cell -> bits))(2) = 0;
    }
    else if (l == 2) {
      (*(cell -> bits))(1) = 0;
      (*(cell -> bits))(2) = 1;
    }
    else if (l == 3) {
      (*(cell -> bits))(1) = 1;
      (*(cell -> bits))(2) = 0;
    }
    else {
      (*(cell -> bits))(1) = 1;
      (*(cell -> bits))(2) = 1;
    }
  }

  // drift one PCM cell
  void drift_one_cell(Cell* cell, long& time){
    cell -> t += time;
    cell -> R = get_drifted_resistance(cell);
    int newlevel = compute_new_level(cell -> R);
    if (newlevel != cell -> L){
      cell -> L = newlevel;
      update_bits_in_cell(cell);
    }
  }

  // drift all the PCM cells
  void drift_all_cells(vector< Cell* >& cells, long& time){
    for(long i = 0; i < (long) cells.size(); i++){
      drift_one_cell(cells[i], time);
    }
  }

		       
  // compute the cell level from the new bits stored in the cell
  int compute_cell_level_from_bits(const Cell* cell){
    int level = 1;
    int exp;
    for(long i = 1; i <= BITS_PER_CELL; i++){
      if (IsOne((*(cell -> bits))(i))){
	exp = BITS_PER_CELL-i;
	level += (int) pow((double) 2, exp);
      }
    }    
  
    return level;
  }

  // // compute bit error rate given the PCM cells.
  // double compute_bit_error_rate(const vector< Cell* >& cells, const mat_GF2& H){
  //   long errnum = 0;
  //   for(long i = 0; i < (long) cells.size(); i++){
  //     Cell* cel = cells[i];
  //     for(long j = 1; j <= BITS_PER_CELL; j++){
  // 	if ((*(cel -> bits))(j) != (*(cel -> bits0))(j))
  // 	  errnum++;
  //     }
  //   }
  //   return ((double) errnum / H.NumCols());
  // }

  // // compute bit error rate given the PCM cells.
  // double compute_bit_error_rate(const vector< Cell* >& cells, const long& n){
  //   long errnum = 0;
  //   for(long i = 0; i < (long) cells.size(); i++){
  //     Cell* cel = cells[i];
  //     for(long j = 1; j <= BITS_PER_CELL; j++){
  // 	if ((*(cel -> bits))(j) != (*(cel -> bits0))(j))
  // 	  errnum++;
  //     }
  //   }
  //   return ((double) errnum / (double) n);
  // }


  // compute the mu and delta for the Gaussian distribution of the
  // resistance after the cells are drifted t time.
  // FIXME: currently we don't use the cached rd trick as we do in the
  // maple code.
  // compute mu and sigma^2
  inline void compute_mean_sdev2_t(double* md, const int& l, const long& t){
    double lgt = log10(t);
    md[0] = resistance_means[l] + exponent_means[l] * lgt;
    md[1] = pow(resistance_dev, (double) 2.0) + pow(exponent_devs[l],(double) 2.0) * pow(lgt, (double) 2.0);
  }

  // evaluate a probability density function where 
  // mean = md[0], and sigma^2 = md[1]
  inline double eval_pdf(const double* md, const double& r){
    return exp((pow((r - md[0]), (double) 2))/(-2*md[1]))/sqrt(2 * M_PI * md[1]);  
  }

  // FIXME: specifically designed for 4-level cells.
  inline double compute_P(const double& r, const long& t, 
			  const int l1, const int l2){  
    double mdl1[2], mdl2[2];
    compute_mean_sdev2_t(mdl1, l1, t);
    compute_mean_sdev2_t(mdl2, l2, t);
    return ((eval_pdf(mdl1, r) + eval_pdf(mdl2, r)) / ((double) MAX_CELL_LEVEL));
  }

  // FIXME: specifically designed for 4-level cells
  void compute_p1_unroll(double* ps, const double& r, const long& t){
    double cp01, cp11;
    double cp02, cp12;
    cp01 = compute_P(r, t, 1, 2);
    cp11 = compute_P(r, t, 3, 4);
    cp02 = compute_P(r, t, 1, 3);
    cp12 = compute_P(r, t, 2, 4);
    ps[0] = cp11/(cp01 + cp11);
    ps[1] = cp12/(cp02 + cp12);
  }

  // compute the initial probabilities for one cell
  void compute_initial_prob(double* ps, const Cell* cel){
    long t = cel -> t;
    int j;

    if (t == 0){
      for(j = 0; j < BITS_PER_CELL; j++){
	if (IsOne((*(cel -> bits))(j+1))){
	  ps[j] = 1;
	}
	else{
	  ps[j] = 0;
	}
      }
    }
    else{
      compute_p1_unroll(ps, cel -> R, t);
    }
  }

  // compute the initial probabiliteis for all the variable nodes/cells
  void compute_initial_probs(double* probs, const vector< Cell* >& cells){
    int index = 1;
    for(long i = 0; i < (long) cells.size(); i++){
      Cell* cel = cells[i];
      double ps[BITS_PER_CELL];
      compute_initial_prob(ps, cel);
      for(int j = 0; j < BITS_PER_CELL; j++, index++){
	probs[index] = ps[j];
      }
    }
    probs[0] = 0.0;
  }

  // compute the initial error probabilities then invoke the LDPC decoding.
  bool error_correct(vec_GF2& correctedBits, int& iternum, const vector< Cell* >& cells, const mat_GF2& H){
    long n = H.NumCols();
    double* initialProbs = new double[n+1];
    compute_initial_probs(initialProbs, cells);
    if (DEBUG){
      cout << "The initial probabilities" << endl;
      for(long i = 0; i < n; i++){
	cout << initialProbs[i] << " ";
      }
    }
    bool convergeflag =  LDPCDecoder::decodeSPA(correctedBits, H, iternum, initialProbs);  
    delete [] initialProbs;
    return convergeflag;
  }

  // check whether the i-th cell contains error
  bool is_faulty_cell(const long& i, vec_GF2& correctedBits, 
		      const vector< Cell* >& cells){
    long bindex;
    for(int j = 1; j <= BITS_PER_CELL; j++){
      bindex = i * BITS_PER_CELL + j;
      if ((*(cells[i] -> bits))(j) != correctedBits(bindex))
	return true;
    }
    return false;
  }

  // updates the faulty bits found by the ECC store the correct bits of
  // the faulty cells into the error log
  void update_error_log(std::set< pair<GF2, long>, ltpair >& errorLog, 
			const vector< Cell* >& cells, 
			vec_GF2& correctedBits) {
    long bindex;
    for(long i = 0; i < NUM_CELL; i++) {
      is_faulty_cell(i, correctedBits, cells);
      if (is_faulty_cell(i, correctedBits, cells)) {
	for(int j = 1; j <= BITS_PER_CELL; j++){
	  bindex = i * BITS_PER_CELL + j;
	  pair<GF2, long> log (correctedBits(bindex), bindex - 1);
	  errorLog.insert(log);
	}
      }
    }
    //   for(int j = 1; j <= BITS_PER_CELL; j++){
    //     bindex = i * BITS_PER_CELL + j;
    //     if ((*(cells[i] -> bits))(j) != correctedBits(bindex)) {
    // 	// update the error log
    // 	//pair<GF2, long> log ((*(cells[i] -> bits))(j), bindex - 1);
    // 	// we memoize correct values of the errors, not errors any more (see above line)
    // 	pair<GF2, long> log (correctedBits(bindex), bindex - 1);
    // 	errorLog.insert(log);
    //     }
    //   }
    // }
  }

  // print out the list of cells.
  void print_cells_bits_only(const vector< Cell* >& cells){
    for(long j = 0; j < (long) cells.size(); j++){
      cout << *(cells[j]->bits0) << " ";
    }
    cout << endl;
  }

  void print_cells(const vector< Cell* >& cells){
    for(long j = 0; j < (long) cells.size(); j++){
      cout << *(cells[j]) << endl;
    }
  }


  // refresh a cell according to the newly stored bits.
  inline void refresh_one_cell(Cell* cel){
    cel -> L = compute_cell_level_from_bits(cel);
    cel -> t = 0;
    cel -> R0 = get_new_resistance(cel -> L);
    cel -> R = cel -> R0;
    cel -> V = get_new_drift_exponent(cel -> L);
  }


  // cleans up all the cells.
  void clean_up_all_cells(const int& i, vector< Cell* >& cells, 
			  const vec_GF2& correctedBits, 
			  std::set< pair<GF2, long>, ltpair >& errorLog,
			  bool remove) {

    for(long j = 1; j <= CELL_PER_SEG; j++) {
      // cindex starts from 0
      long cindex = (i - 1) * CELL_PER_SEG + j - 1;
      for(long k = 1; k <= BITS_PER_CELL; k++) {
	// FIXME: check here cell 0 -> bits 1, 2
	long bindex = cindex * BITS_PER_CELL + k;
	// correct the faulty bits.
	(*(cells[cindex] -> bits))(k) = correctedBits(bindex);
	// remove the entry from errorLog if this bit is found
	// faulty by ECC
	if (use_error_log_before_ecc) {
	  if (remove){
	    std::set< pair<GF2, long>, ltpair >::iterator it;
	    pair<GF2, long> key((*(cells[cindex] -> bits))(k), bindex - 1);
	    it = errorLog.find(key);
	    if (it != errorLog.end()){
	      errorLog.erase(it);
	    }
	  }
	}
      }
      // refresh this cell
      // update the resistance, exponent from the newly computed level.
      // The reset elapsed time to 0.
      refresh_one_cell(cells[cindex]);
    }
  }


  // cleans up all the cells 
  // used in the naive scheme
  void clean_up_all_cells_basic(vector< Cell* >& cells, const vec_GF2& correctedBits){
    for(long cindex = 0; cindex < (long) cells.size(); cindex++){
      for(long k = 1; k <= BITS_PER_CELL; k++) {
	// FIXME: check here cell 0 -> bits 1, 2
	long bindex = cindex * BITS_PER_CELL + k;
	// correct the faulty bits.
	(*(cells[cindex] -> bits))(k) = correctedBits(bindex);
      }
      // refresh this cell
      // update the resistance, exponent from the newly computed level.
      // The reset elapsed time to 0.
      refresh_one_cell(cells[cindex]);
    }
  }

  // cleans up the cells if their stored bits are faulty as recorded in
  // the error log.
  void correct_if_in_log(const std::set< pair<GF2, long>, ltpair >& errorLog,
			 int& correctedBeforeECC, 
			 const long& k, const long& bindex,
			 Cell* cell, 
			 vector< Cell* >& to_refresh){  
    std::set< pair<GF2, long>, ltpair>::iterator iter;
    for(iter = errorLog.begin(); iter != errorLog.end(); iter++){
      if ((*iter).second == bindex){
	// if the bits do not equal to the correct bits
	if ((*iter).first != (*(cell -> bits))(k)){
	  correctedBeforeECC++;
	  (*(cell -> bits))(k) += 1;
	  to_refresh.push_back(cell);
	}
      }
    }
  }

  // cleans up the segment i of the cells using the error log.
  void cleans_up_segment_with_error_log(const int& i, vector< Cell* >& cells,  
					std::set< pair<GF2, long>, ltpair >& errorLog, 
					vector< Cell* >& cells_to_refresh, 
					int& correctedBeforeECC){
    if (errorLog.size() > 0){
      for(long j = 1; j <= CELL_PER_SEG; j++) {
	// cindex starts from 0
	long cindex = (i - 1) * CELL_PER_SEG + j - 1;
	for(long k = 1; k <= BITS_PER_CELL; k++) {
	  //FIXME: check here cell 0 -> bits 0, 1
	  long bindex = cindex * BITS_PER_CELL + k - 1;
	  correct_if_in_log(errorLog, correctedBeforeECC,  
			    k, bindex, cells[cindex], cells_to_refresh);
	}
      }
    }
  }


  // update performance result obtained from one round of experiment
  // specific for LDPC
  void update_common_global_performance_counters(const long& totalIters, const long& nofaultybits, 
						 const long& correctblkno, const long& totalErrorBlocks, 
						 const long& undetblkno, const long& detblkno){
    average_iter_num += totalIters;
    average_error_num += nofaultybits;
    average_correct_blk_num += correctblkno;
    average_error_blk_num += totalErrorBlocks;
    average_undet_blk_num += undetblkno;
    average_det_blk_num += detblkno;
  }

  // fill vector v with x
  void fill_vec_GF2(vec_GF2& v, int x){
    long n = v.length();
    for(long i = 1; i <= n; ++i){
      v(i) = to_GF2(x);
    }
  }

  // pad the metadata with random bits so that metadata is long enough
  // to be encoded by the LDPC code C2 (N2, K2)
  void pad_metadata(vec_GF2& metadata, long& numpadbits){
    if (metadata.length() > K2){
      cout << "Error: the unpad metadata is too long for the current choose of  H2!" << endl;
      exit(EXIT_FAILURE);
    }
    numpadbits = K2 - metadata.length();
    // seed NTL's random vector generator.
    srand(time(NULL) * rand());
    SetSeed(to_ZZ(rand()));
    append(metadata, random_vec_GF2(K2 - metadata.length()));
  }

  // store the error log into an another array of cells.  We assume that
  // the Huffman lookup table should be stored in single level cells so
  // that it won't be effected by cell level drifting.
  void save_error_log_in_mcells(vector< Cell* >& mcells, 
				const std::set< pair<GF2, long>, ltpair >& errorLog, 
				std::vector<TreeNode*>& trees, 
				const mat_GF2& H2,       
				const mat_GF2& G2,
				vec_GF2& metadata,        // stores the compressed version of metabits
				vec_GF2& encodedMetadata, // the encoded version of the above metadata 
				long& numpadbits){
    encodedMetadata.kill();
    metadata.kill();
    // cleans up the previous data stored in mcells.
    if (mcells.size() != 0){
      free_PCM_cells(mcells);
      mcells.clear();
    }
    vec_GF2 metabits;            // a bit vector holding the correct values of the faulty bits from errorLog 
    vector<int> levels;          
    //bool appendToStore = false;  // true if the metadata are further appended with one more bit

    // 1. update metabits with the correct values of errors
    // fill other slots in metabits with 1, this is to indicator the relevance of those 1's
    metabits.SetLength(N);
    fill_vec_GF2(metabits, 1);
    std::set< pair<GF2, long>, ltpair >::iterator it;
    for(it = errorLog.begin(); it != errorLog.end(); it++){
      metabits((*it).second + 1) = (*it).first;
    } 

    // 2. compute the cell levels represented by metabits
    CellUtil::compute_cell_levels(levels, metabits);

    // 3. compress the levels represented by metabits using Huffman code
    Huffman::huffman_encode(metadata, levels, trees);
#ifdef VERBOSE
    cout << "Size of unpaded metadata:" << metadata.length() << endl;
    cout << "Compression rate:" << (double) metadata.length() / (double) N << endl;
#endif

    // update performance counter
    average_compressed_metadata_size += metadata.length();

    // pad the metadata with random bits so that metadata.length() == K2
    // the random padding is to increase the performance of LDPC soft decoding
    // i.e. for a more precise estimation of initial probabilities.
    pad_metadata(metadata, numpadbits);

    // encode metadata using LDPC code C2 (N2, K2)
    encodedMetadata.SetLength(N2);
    LDPCEncoder::encode(encodedMetadata, metadata, G2);
#ifdef VERBOSE
    cout << "num pad bits:" << numpadbits << endl;
    cout << "Size of encoded metadata:" << encodedMetadata.length() << endl;
    //cout << encodedMetadata << endl;
#endif

    // append one more bit if the length of metadata is not even.
    //   if ((encodedMetadata.length() % 2) != 0){
    //     append(encodedMetadata, encodedMetadata(encodedMetadata.length()));
    //     appendToStore = true;
    // #ifdef VERBOSE
    //     cout << "The encoded metadata has been padded with one more bit" << endl;
    // #endif
    //   }

    // store the compressed data into PCM cells
    initialize_PCM_cells(mcells, encodedMetadata);
  }


  // compute the corresponding bits from a given cell level
  void compute_bit_from_level(vec_GF2& cellbits, const int& l){
    if (l == 1) {
      cellbits(1) = 0;
      cellbits(2) = 0;
    }
    else if (l == 2) {
      cellbits(1) = 0;
      cellbits(2) = 1;
    }
    else if (l == 3) {
      cellbits(1) = 1;
      cellbits(2) = 0;
    }
    else {
      cellbits(1) = 1;
      cellbits(2) = 1;
    }
  }

  // compute the corresponding bits from the levels of the cells
  void compute_bits_from_levels(vec_GF2& metabits, const vector<int>& levels){
    vec_GF2 cellbits;
    cellbits.SetLength(BITS_PER_CELL);
    int level;
    for(long i = 0; i < (long) levels.size(); ++i){
      level = levels[i];
      compute_bit_from_level(cellbits, level);
      append(metabits, cellbits);
    }
  }

  // checks whether the bits in a cell are all ones
  bool contain_all_ones(long& i, const vec_GF2& metabits){
    for(long j = 0; j < BITS_PER_CELL; ++j){
      if (IsZero(metabits(i+j) ))
	return false;
    }
    return true;
  }


  // create error log from metabits
  void create_error_log_from_metabits(std::set< pair<GF2, long>, ltpair >& errorLog,
				      const vec_GF2& metabits){
    long bindex;
    for(long i = 1; i <= metabits.length(); i+=BITS_PER_CELL){
      if (!contain_all_ones(i, metabits)){
	for(long j = 0; j < BITS_PER_CELL; ++j){
	  bindex = i + j;
	  pair<GF2, long> log (metabits(bindex), bindex - 1);
	  errorLog.insert(log);
	}
      }
    }
  }


  // update performance counters on the metadata decoding
  void update_performance_counters_of_metadata_decoder(const bool& successflag, 
						       const vector< Cell* >& mcells,
						       const vec_GF2& correctedBits,
						       const vec_GF2& encodedMetadata,
						       const int& iternum, 
						       const mat_GF2& H2){
    long numerrs = 0;
    trial_num_meta++;
    average_ber_before_decode_meta += CellUtil::compute_bit_error_rate(mcells, H2.NumCols());
    average_iter_num_meta += iternum;
    numerrs = weight(encodedMetadata - correctedBits);
    average_error_num_meta += numerrs;
    if (numerrs == 0)
      average_correct_blk_num_meta++;
    else
      average_error_blk_num_meta++;
  }

  void free_code_tree(TreeNode* tree){
    if (tree -> leftChild){
      free_code_tree(tree -> leftChild);
    }
    if (tree -> rightChild){
      free_code_tree(tree -> rightChild);
    }
    delete tree;
  }

  void free_trees(std::vector<TreeNode*>& trees){
    if (trees.size() > 0){
      free_code_tree(trees[0]);
      trees.clear();
    }
  }

  // load from mcells to create error log
  // this is a procedure opposite to save_error_log_in_mcells
  void retrieve_error_log_from_mcells(std::set< pair<GF2, long>, ltpair >& errorLog,
				      const vector< Cell* >& mcells, 
				      const long& numpadbits, 
				      std::vector<TreeNode*>& trees,
				      vec_GF2& metadata,
				      vec_GF2& encodedMetadata,
				      const mat_GF2& H2){
    vec_GF2 correctedBits;
    vector<int> levels;
    vec_GF2 metabits;
    bool successflag;
    if (errorLog.size() > 0){
      vec_GF2 correctedBits;
      int iternum = 0;
      errorLog.clear();
#ifdef VERBOSE
      cout << "Correcting the metadata" << endl;
#endif

      successflag = error_correct(correctedBits, iternum, mcells, H2);
      update_performance_counters_of_metadata_decoder(successflag, mcells, 
						      correctedBits, encodedMetadata, 
						      iternum, H2);
      if (successflag) {
	// unpad the metadata
	metadata.SetLength(K2 - numpadbits);
	// decompress using Huffman decoding
	Huffman::huffman_decode(levels, metadata, trees);
	// free the tree node
	free_trees(trees);

	// compute metabits 
	compute_bits_from_levels(metabits, levels);
	// create error log from metabits
	create_error_log_from_metabits(errorLog, metabits);
      }
    }
  }

#ifndef PR1   // if proposal 1  is not used

  // periodically scrub the memory cells
  // codeword: a randomly generated codeword to be stored in the cells.
  // H: parity check matrix for LDPC code
  void multiphase_scrubbing_ldpc(const vec_GF2& codeword, const mat_GF2& H, 
			    const mat_GF2& H2,       const mat_GF2& G2){
    vector< pair< double, double> > bers;   
    vector< Cell* > cells;                  // cells for storing payload data. index starts from 0
    vector< Cell* > mcells;                 // cells for storing the meta data. index starts from 0
    std::set< pair<GF2, long>, ltpair > errorLog;   // the error log
    std::vector<TreeNode*> trees;          // used for storing huffman tree for compressing meta data
    vec_GF2 metadata; 
    vec_GF2 encodedMetadata;                // the encoded version of the above metadata

    long t;                                 // current time t
    int i = 1;                              // index of the current segment being scrubbed.
    int s = NUM_CELL/CELL_PER_SEG;          // number of segments which the cells are divided into.
    long gap = (long) scrub_gap;
  
    double berbefore;                       // BER before applying EC
    double berafter;                        // BER after applying EC
    bool successflag;                       // indicates whether a decoding converges
    long nonconvergenum = 0;                // counting consecutive nonconverge decodings
    long bitsRefreshedBeforeEC = 0;         // counting the number of bits refreshed by error log
    long bitsMemoized = 0;                  // counting the size of error log.
    long numpadbits = 0;                    // the number of bits which are padded

    // initializing performance counters
    long correctblkno     = 0;
    long undetblkno       = 0;
    long detblkno         = 0;
    long totalErrorBlocks = 0;
    long totalIters       = 0;
    long nosegscrubbed    = 0;
    long nofaultybits     = 0;

    // initialization of cells
    initialize_PCM_cells(cells, codeword);

    if (DEBUG){
      cout << "The initialized cells:" << endl;
      print_cells(cells);
    }
  
    t = gap;
    drift_all_cells(cells, t);
  
    if (DEBUG){
      cout << "The drifted cells" << endl;
      print_cells(cells);
    }
  
#ifdef VERBOSE
    if (use_error_log_before_ecc) 
      cout << "Use error log before EC" << endl;
    else
      cout << "Don't use error log" << endl;
    cout << "Total number of segments: " << s << endl;
#endif

    while (t <= total_time) {

#ifdef VERBOSE
      cout << endl << "Segment#: " << i << endl;
      cout << "Time: " << t << endl;
#endif

      if (use_error_log_before_ecc) {
	vector< Cell* > cells_to_refresh;
	int correctedBeforeECC = 0;

	retrieve_error_log_from_mcells(errorLog, mcells, numpadbits, trees, metadata, encodedMetadata, H2);
	cleans_up_segment_with_error_log(i, cells, errorLog, cells_to_refresh, correctedBeforeECC);

	// update two local counters
	bitsRefreshedBeforeEC += correctedBeforeECC; 
	bitsMemoized += errorLog.size();

#ifdef VERBOSE
	cout << "Number of faulty bits corrected before EC: " << correctedBeforeECC << endl;
	cout << "Total number of bits in error log: " << errorLog.size() << endl;
#endif

	//std::unique(cells_to_refresh.begin(), cells_to_refresh.end());
	for(long j = 0; j < (long) cells_to_refresh.size(); j++){
	  refresh_one_cell(cells_to_refresh[j]);
	}
      
#ifdef VERBOSE
	cout << "Number of cells to refresh:" << cells_to_refresh.size() << endl;
#endif
      }
      //    else {

      //#ifdef VERBOSE
      //cout << "Decoding result is not reliable, while error log still has correct bits, we don't do anything after this scrubbing" << endl;
      //#endif
      //    }

      // save the ber before ECC.
      berbefore = CellUtil::compute_bit_error_rate(cells, H.NumCols());
#ifdef VERBOSE
      cout << "BER before EC:" << berbefore << endl;   
      cout << "Correcting payload data" << endl;
#endif

      // step 2. invoke ECC (LDPC here)
      vec_GF2 correctedBits;
      int iternum = 0;
      successflag = error_correct(correctedBits, iternum, cells, H);

      if (successflag) 
	nonconvergenum = 0;
      else
	nonconvergenum++;

      // if the decoding is successful we refresh
      // or if the decodig is not terminating, and the error log is
      // empty i.e., there's no more good hope, we use best guess.
      if (successflag || (errorLog.size() == 0)){
	if (use_error_log_before_ecc) {
	  errorLog.clear();
	  update_error_log(errorLog, cells, correctedBits);
	}

	// clean up all the cells in the current
	// segment by writing the bits returned by LDPC decoding
	// back to the segment      
	clean_up_all_cells(i, cells, correctedBits, errorLog, true);
      }
      else if ((!successflag) && (nonconvergenum >= (s-1))){
#ifdef VERBOSE
	cout << "Consecutive nonconverge decoding times: " << nonconvergenum << endl;
	cout << "Using best guess to clean up" << endl;
#endif
	clean_up_all_cells(i, cells, correctedBits, errorLog, false);
      }
      else{
#ifdef VERBOSE
	cout << "Nonconverge, do nothing!" << endl;
#endif
      }

      // write the error log into mcells
      save_error_log_in_mcells(mcells, errorLog, trees, H2, G2, metadata, encodedMetadata, numpadbits);

      // save the ber after EC
      berafter = CellUtil::compute_bit_error_rate(cells, H.NumCols());
#ifdef VERBOSE
      cout << "BER after EC: " << berafter << endl;
#endif

      // update the global bit error rate vector
      bit_error_rates[nosegscrubbed * 2] += berbefore;
      bit_error_rates[nosegscrubbed * 2 + 1] += berafter;

      totalIters += iternum;
      nosegscrubbed++;
      nofaultybits = weight(codeword - correctedBits);
      if (nofaultybits == 0) {
	if (successflag)
	  correctblkno++;
      }
      else {
	totalErrorBlocks++;
	if (successflag)
	  undetblkno++;
	else
	  detblkno++;
      }

      // increase time and move to the next segment (by increasing
      // the index i)
      // if this is the last segment, go back to scrub the first one   
      if (i == s) 
	i = 1;
      else
	i++;
      t+= gap;
      // drift all the cells containing payload data
      drift_all_cells(cells, gap);
      // drift all the cells containing metadata
      drift_all_cells(mcells, gap);
    }
  
    // update performance counter
    num_seg_scrubbed += nosegscrubbed;
    average_cleanup_bit_num += bitsRefreshedBeforeEC; 
    average_mem_bit_num += bitsMemoized;
    update_common_global_performance_counters(totalIters, nofaultybits, correctblkno, 
					      totalErrorBlocks, undetblkno, detblkno);
    
    free_PCM_cells(cells);
    free_PCM_cells(mcells);
    free_trees(trees);
  }

#else
  // updates the correct bits of the next segment found by the ECC
  void update_error_log_with_next_segment(vec_GF2& errorLog, 
					  vec_GF2& correctedBits, 
					  int& nextseg) {
    long bindex;
    long k = 1;
    for(long i = (nextseg - 1) * CELL_PER_SEG; i < nextseg * CELL_PER_SEG; i++){
      for(int j = 1; j <= BITS_PER_CELL; j++){
	bindex = i * BITS_PER_CELL + j;
	errorLog(k) = correctedBits(bindex);
	k++;
      }
    }
  }


  // count the number of errors found in segment seg
  long count_error_num_in_seg(vector< Cell* >& cells,
			      const int& seg, 
			      const vec_GF2& correctedBits) {
    long result = 0;
    long bindex;
    for(long i = (seg - 1) * CELL_PER_SEG; i < seg * CELL_PER_SEG; i++){
      for(int j = 1; j <= BITS_PER_CELL; j++){
	bindex = i * BITS_PER_CELL + j;
	if ((*(cells[i] -> bits))(j) != correctedBits(bindex)) {
	  result++;
	}
      }
    }
    return result;
  }

  // find the most errorneous segment excluding nextseg.
  int find_most_errorneous_segment(vector< Cell* >& cells,
				   const int& nextseg, 
				   const vec_GF2& correctedBits) {
    int result = 0;
    long maxerrornum = 0;
    long errornum;
    for(int i = 1; i <= (int) (NUM_CELL/CELL_PER_SEG); i++){
      if (i != nextseg){
	errornum = count_error_num_in_seg(cells, i, correctedBits);
	if (errornum > maxerrornum){
	  maxerrornum = errornum;
	  result = i;
	}
      }
    }
#ifdef VERBOSE
    cout << "The most errorneous segment:" << result;
    cout << " Error#:" << maxerrornum << endl;
#endif
    return result;
  }


  // cleans up the most erroneous segment of cells.
  void clean_up_most_erroneous_segment(vector< Cell* >& cells,
				       const int& nextseg, 
				       const vec_GF2& correctedBits) {
    int i = find_most_errorneous_segment(cells, nextseg, correctedBits);
    long cindex;
    for(long j = 1; j <= CELL_PER_SEG; j++) {
      // cindex starts from 0
      cindex = (i - 1) * CELL_PER_SEG + j - 1;
      for(long k = 1; k <= BITS_PER_CELL; k++) {
	// FIXME: check here cell 0 -> bits 1, 2
	long bindex = cindex * BITS_PER_CELL + k;
	// correct the faulty bits.
	(*(cells[cindex] -> bits))(k) = correctedBits(bindex);
      }
      // refresh this cell
      // update the resistance, exponent from the newly computed level.
      // The reset elapsed time to 0.
      refresh_one_cell(cells[cindex]);
    }
  }


  // cleans up the segment i of the cells using the error log.
  void cleans_up_segment_with_error_log_prop1(const int& i, vector< Cell* >& cells,  
					      const vec_GF2& errorLog, 
					      vector< Cell* >& cells_to_refresh, 
					      int& correctedBeforeECC){
    long bindex = 1;
    if (!IsZero(errorLog)){
      for(long j = 1; j <= CELL_PER_SEG; j++) {
	// cindex starts from 0
	long cindex = (i - 1) * CELL_PER_SEG + j - 1;
	for(long k = 1; k <= BITS_PER_CELL; k++) {
	  //FIXME: check here cell 0 -> bits 1, 2
	  (*(cells[cindex] -> bits))(k) = errorLog(bindex);
	  cells_to_refresh.push_back(cells[cindex]);
	  correctedBeforeECC++;
	  bindex++;
	}
      }
    }
  }


  // periodically scrub the memory cells 
  // codeword: a randomly generated codeword to be stored in the cells.
  // H: parity check matrix for LDPC code
  void multiphase_scrubbing_ldpc(const vec_GF2& codeword, const mat_GF2& H){
    vector< pair< double, double> > bers;   
    vector< Cell* > cells;                  // cells for storing payload data. index starts from 0

    vec_GF2 errorLog;
    errorLog.SetLength(BITS_PER_CELL * CELL_PER_SEG);

    long t;                                 // current time t
    int i = 1;                              // index of the current segment being scrubbed.
    int nextseg;
    int s = NUM_CELL/CELL_PER_SEG;          // number of segments which the cells are divided into.
    long gap = (long) scrub_gap;

    double berbefore;                       // BER before applying EC
    double berafter;                        // BER after applying EC
    bool successflag;                       // indicates whether a decoding converges
    long nonconvergenum = 0;                // counting consecutive nonconverge decodings
    long bitsRefreshedBeforeEC = 0;         // counting the number of bits refreshed by error log
    long bitsMemoized = 0;                  // counting the size of error log.

    // initializing performance counters
    long correctblkno     = 0;
    long undetblkno       = 0;
    long detblkno         = 0;
    long totalErrorBlocks = 0;
    long totalIters       = 0;
    long nosegscrubbed    = 0;
    long nofaultybits     = 0;

    // initialization of cells
    initialize_PCM_cells(cells, codeword);

    if (DEBUG){
      cout << "The initialized cells:" << endl;
      print_cells(cells);
    }
  
    t = gap;
    drift_all_cells(cells, t);
  
    if (DEBUG){
      cout << "The drifted cells" << endl;
      print_cells(cells);
    }
  
#ifdef VERBOSE
    if (use_error_log_before_ecc) 
      cout << "Use error log before EC" << endl;
    else
      cout << "Don't use error log" << endl;
    cout << "Total number of segments: " << s << endl;
#endif

    while (t <= total_time) {

#ifdef VERBOSE
      cout << endl << "Segment#: " << i << endl;
      cout << "Time: " << t << endl;
#endif

      nextseg = i + 1;
      if (i == s) {
	nextseg = 1;
      }

      if (use_error_log_before_ecc) {
	vector< Cell* > cells_to_refresh;
	int correctedBeforeECC = 0;
	// FIXME: need to modify the following function. 
	// the preclean-up uses the *correct* bits stored in the error log to refresh the whole segment.
	cleans_up_segment_with_error_log_prop1(i, cells, errorLog, cells_to_refresh, correctedBeforeECC);

	// update two local counters
	bitsRefreshedBeforeEC += correctedBeforeECC; 
	bitsMemoized += CELL_PER_SEG * BITS_PER_CELL;

#ifdef VERBOSE
	cout << "Number of faulty bits corrected before EC: " << correctedBeforeECC << endl;
	cout << "Total number of bits in error log: " << CELL_PER_SEG * BITS_PER_CELL << endl;
#endif
	for(long j = 0; j < (long) cells_to_refresh.size(); j++){
	  refresh_one_cell(cells_to_refresh[j]);
	}
#ifdef VERBOSE
	cout << "Number of cells to refresh:" << cells_to_refresh.size() << endl;
#endif
      }
      else {
#ifdef VERBOSE
	cout << "Decoding result is not reliable, while error log still has correct bits, we don't do anything after this scrubbing" << endl;
#endif
      }

      // save the ber before ECC.
      berbefore = CellUtil::compute_bit_error_rate(cells, H.NumCols());
#ifdef VERBOSE
      cout << "BER before EC:" << berbefore << endl;   
#endif

      // step 2. invoke ECC (LDPC here)
      vec_GF2 correctedBits;
      int iternum = 0;
      successflag = error_correct(correctedBits, iternum, cells, H);

      if (successflag) 
	nonconvergenum = 0;
      else
	nonconvergenum++;

      // we refresh the segment which is most errorneous excluding the
      // next segment.
#ifdef VERBOSE
      if (!successflag)
	cout << "Nonconverge, using best guess to clean up and update error log" << endl;
#endif
      if (use_error_log_before_ecc) {
	update_error_log_with_next_segment(errorLog, correctedBits, nextseg);
      }
      clean_up_most_erroneous_segment(cells, nextseg, correctedBits);

      // save the ber after EC
      berafter = CellUtil::compute_bit_error_rate(cells, H.NumCols());
#ifdef VERBOSE
      cout << "BER after EC: " << berafter << endl;
#endif

      // update the global bit error rate vector
      bit_error_rates[nosegscrubbed * 2] += berbefore;
      bit_error_rates[nosegscrubbed * 2 + 1] += berafter;

      totalIters += iternum;
      nosegscrubbed++;
      nofaultybits = weight(codeword - correctedBits);
      if (nofaultybits == 0) {
	if (successflag)
	  correctblkno++;
      }
      else {
	totalErrorBlocks++;
	if (successflag)
	  undetblkno++;
	else
	  detblkno++;
      }
      // increase time and move to the next segment (by increasing
      // the index i)
      // if this is the last segment, go back to scrub the first one   
      i = nextseg;
      t+= gap;
      // drift all the cells.
      drift_all_cells(cells, gap);
    }
  
    // update performance counter
    num_seg_scrubbed += nosegscrubbed;
    average_cleanup_bit_num += bitsRefreshedBeforeEC; 
    average_mem_bit_num += bitsMemoized;
    update_common_global_performance_counters(totalIters, nofaultybits, correctblkno, 
					      totalErrorBlocks, undetblkno, detblkno);
    free_PCM_cells(cells);
  }
#endif // PR1

  // simulating the scheme without any scrubbing
  void static_sensing_ldpc(const vec_GF2& codeword, const mat_GF2& H){
    vector< Cell* > cells; 
    double initialProbs[N];
    vec_GF2 correctedBits;
    int iternum = 0;
    long t = total_time;
    initialize_PCM_cells(cells, codeword);
    drift_all_cells(cells, t);
    compute_initial_probs(initialProbs, cells);
#ifdef VERBOSE
    cout << "BER before EC: " << CellUtil::compute_bit_error_rate(cells, H.NumCols()) << endl;
#endif
    LDPCDecoder::decodeSPA(correctedBits, H, iternum, initialProbs);
    free_PCM_cells(cells);
  }

		  
  void naive_scrubbing_ldpc(const vec_GF2& codeword, const mat_GF2& H){
    vector< pair< double, double> > bers;   
    vector< Cell* > cells;                  // index starts from 0
    long t;                                 // current time t
    long gap = scrub_gap_naive;

    double berbefore; // BER before applying EC
    double berafter;  // BER after applying EC
    bool successflag; // indicates whether a decoding converges

    // initializing performance counters
    long correctblkno     = 0;
    long undetblkno       = 0;
    long detblkno         = 0;
    long totalErrorBlocks = 0;
    long totalIters       = 0;
    long noscrubbed    = 0;
    long nofaultybits     = 0;

    // initialization of cells
    initialize_PCM_cells(cells, codeword);

    if (DEBUG){
      cout << "The initialized cells:" << endl;
      print_cells(cells);
    }
  
    t = gap;
    drift_all_cells(cells, t);
  
    if (DEBUG){
      cout << "drifted cells" << endl;
      print_cells(cells);
    }
  
    while (t <= total_time) {
#ifdef VERBOSE
      cout << "Time: " << t << endl;
#endif

      // save the ber before ECC.
      berbefore = CellUtil::compute_bit_error_rate(cells, H.NumCols());
#ifdef VERBOSE
      cout << "BER before EC:" << berbefore << endl;   
#endif

      // step 2. invoke ECC (LDPC here)
      vec_GF2 correctedBits;
      int iternum = 0;
      successflag = error_correct(correctedBits, iternum, cells, H);

      // clean up all the cells in the current
      // segment by writing the bits returned by LDPC decoding
      // back to the segment      
      clean_up_all_cells_basic(cells, correctedBits);

      // save the ber after ECC.
      berafter = CellUtil::compute_bit_error_rate(cells, H.NumCols());
#ifdef VERBOSE
      cout << "BER after EC: " << berafter << endl;
#endif

      // update the global bit error rate vector
      bit_error_rates[noscrubbed * 2] += berbefore;
      bit_error_rates[noscrubbed * 2 + 1] += berafter;

      totalIters += iternum;
      noscrubbed++;
      nofaultybits = weight(codeword - correctedBits);
      if (nofaultybits == 0) {
	if (successflag)
	  correctblkno++;
      }
      else {
	totalErrorBlocks++;
	if (successflag)
	  undetblkno++;
	else
	  detblkno++;
      }

      // increase time
      t+= gap;
      // drift all the cells.
      drift_all_cells(cells, gap);
    }
  
    // update performance counter
    update_common_global_performance_counters(totalIters, nofaultybits, correctblkno, 
					      totalErrorBlocks, undetblkno, detblkno);
    free_PCM_cells(cells);
  }

  // update performance result obtained from one round of experiment
  // specific for LDPC
  void update_common_global_performance_counters_bch(const long& nofaultybits, 
						     const long& correctblkno, 
						     const long& totalErrorBlocks){
    average_error_num += nofaultybits;
    average_correct_blk_num += correctblkno;
    average_error_blk_num += totalErrorBlocks;
  }



  // naive scrubbing scheme using BCH codes
  void naive_scrubbing_bch(const vec_GF2& codeword){
     vector< pair< double, double> > bers;   
     vector< Cell* > cells;                  // index starts from 0
     long t;                                 // current time t
     long gap = scrub_gap_naive;

     double berbefore; // BER before applying EC
     double berafter;  // BER after applying EC

     // initializing performance counters
     long correctblkno     = 0;
     long totalErrorBlocks = 0;
     long noscrubbed       = 0;
     long nofaultybits     = 0;

     // ugly tricks for BCH's odd block length
     GF2 padBit = codeword(codeword.length());  // the bit used for padding codeword will never go wrong. 
     vec_GF2 codeword_copy = codeword;

     // initialization of cells
     append(codeword_copy, padBit);
     initialize_PCM_cells(cells, codeword_copy);


     if (DEBUG){
       cout << "The initialized cells:" << endl;
       print_cells(cells);
     }

     t = gap;
     drift_all_cells(cells, t);
  
     if (DEBUG){
       cout << "drifted cells" << endl;
       print_cells(cells);
     }
  
     while (t <= total_time) {
#ifdef VERBOSE
       cout << "Time: " << t << endl;
#endif

       // save the ber before ECC.
       berbefore = CellUtil::compute_bit_error_rate(cells, BCHn);
#ifdef VERBOSE
       cout << "BER before EC:" << berbefore << endl;   
#endif

       vec_GF2 correctedBits;
       vec_GF2 noisy_codeword;
       CellUtil::read_bits_from_cell(noisy_codeword, cells);
       if (BCHn * berbefore <= BCHt){
	 // step 2. invoke ECC (BCH hard decoding)
	 // depad noisy_codeword
	 noisy_codeword.SetLength(BCHn);
	 BCHDecoder::decode(correctedBits, noisy_codeword);
	 append(correctedBits, padBit);
	 // clean up all the cells in the current
	 // segment by writing the bits returned by LDPC decoding
	 // back to the segment      
	 clean_up_all_cells_basic(cells, correctedBits);
       }
      else{
	correctedBits = noisy_codeword;
#ifdef VERBOSE
	std::cout << "Exceed BCH codes' Error Correction Capability!" << std::endl;
#endif
      }
       // depad correctedBits
       correctedBits.SetLength(BCHn);

       // save the ber after ECC.
       berafter = CellUtil::compute_bit_error_rate(cells, BCHn);
#ifdef VERBOSE
       cout << "BER after EC: " << berafter << endl;
#endif
       // update the global bit error rate vector
       bit_error_rates[noscrubbed * 2] += berbefore;
       bit_error_rates[noscrubbed * 2 + 1] += berafter;
       noscrubbed++;
       nofaultybits = weight(codeword - correctedBits);
       if (nofaultybits == 0) 
	 correctblkno++;
       else 
	 totalErrorBlocks++;

       // increase time
       t+= gap;
       // drift all the cells.
       drift_all_cells(cells, gap);
     }
  
     // update performance counter
     update_common_global_performance_counters_bch(nofaultybits, correctblkno, totalErrorBlocks);

    free_PCM_cells(cells);
  }


  // multiphase scrubbing with BCH codes.
  void multiphase_scrubbing_bch(const vec_GF2& codeword, const mat_GF2& H2, const mat_GF2& G2){
    vector< pair< double, double> > bers;   
    vector< Cell* > cells;                  // cells for storing payload data. index starts from 0
    vector< Cell* > mcells;                 // cells for storing the meta data. index starts from 0
    std::set< pair<GF2, long>, ltpair > errorLog;   // the error log
    std::vector<TreeNode*> trees;          // used for storing huffman tree for compressing meta data
    vec_GF2 metadata; 
    vec_GF2 encodedMetadata;                // the encoded version of the above metadata

    long t;                                 // current time t
    int i = 1;                              // index of the current segment being scrubbed.
    int s = NUM_CELL/CELL_PER_SEG;          // number of segments which the cells are divided into.
    long gap = (long) scrub_gap;
  
    double berbefore;                       // BER before applying EC
    double berafter;                        // BER after applying EC
    long bitsRefreshedBeforeEC = 0;         // counting the number of bits refreshed by error log
    long bitsMemoized = 0;                  // counting the size of error log.
    long numpadbits = 0;                    // the number of bits which are padded to the metadata


    // initializing performance counters
    long correctblkno     = 0;
    long totalErrorBlocks = 0;
    long nosegscrubbed    = 0;
    long nofaultybits     = 0;  


    // ugly tricks for BCH's odd block length
    GF2 padBit = codeword(codeword.length());  // the bit used for padding codeword will never go wrong. 
    vec_GF2 codeword_copy = codeword;

    // initialization of cells
    append(codeword_copy, padBit);
    initialize_PCM_cells(cells, codeword_copy);

    if (DEBUG){
      cout << "The initialized cells:" << endl;
      print_cells(cells);
    }
  
    t = gap;
    drift_all_cells(cells, t);
  
    if (DEBUG){
      cout << "The drifted cells" << endl;
      print_cells(cells);
    }
  
#ifdef VERBOSE
    if (use_error_log_before_ecc) 
      cout << "Use error log before EC" << endl;
    else
      cout << "Don't use error log" << endl;
    cout << "Total number of segments: " << s << endl;
#endif    


    while (t <= total_time){

#ifdef VERBOSE
      cout << endl << "Segment#: " << i << endl;
      cout << "Time: " << t << endl;
#endif

      if (use_error_log_before_ecc) {
	vector< Cell* > cells_to_refresh;
	int correctedBeforeECC = 0;

	retrieve_error_log_from_mcells(errorLog, mcells, numpadbits, trees, metadata, encodedMetadata, H2);
	cleans_up_segment_with_error_log(i, cells, errorLog, cells_to_refresh, correctedBeforeECC);

	// update two local counters
	bitsRefreshedBeforeEC += correctedBeforeECC; 
	bitsMemoized += errorLog.size();

#ifdef VERBOSE
	cout << "Number of faulty bits corrected before EC: " << correctedBeforeECC << endl;
	cout << "Total number of bits in error log: " << errorLog.size() << endl;
#endif
	for(long j = 0; j < (long) cells_to_refresh.size(); j++){
	  refresh_one_cell(cells_to_refresh[j]);
	}	
#ifdef VERBOSE
	cout << "Number of cells to refresh:" << cells_to_refresh.size() << endl;
#endif
    }

      // save the ber before ECC.
      berbefore = CellUtil::compute_bit_error_rate(cells, BCHn);
#ifdef VERBOSE
      cout << "BER before EC:" << berbefore << endl;   
      cout << "Correcting payload data" << endl;
#endif

      vec_GF2 correctedBits;
      vec_GF2 noisy_codeword;
      CellUtil::read_bits_from_cell(noisy_codeword, cells);
      if (BCHn * berbefore <= BCHt) {
	// step 2. invoke ECC (BCH hard decoding)
	// depad noisy_codeword
	noisy_codeword.SetLength(BCHn);
	BCHDecoder::decode(correctedBits, noisy_codeword);
	append(correctedBits, padBit);

     
	if (use_error_log_before_ecc) {
	  errorLog.clear();
	  update_error_log(errorLog, cells, correctedBits);
	}

	// clean up all the cells in the current
	// segment by writing the bits returned by LDPC decoding
	// back to the segment      
	clean_up_all_cells(i, cells, correctedBits, errorLog, true);

	// depad correctedBits
	correctedBits.SetLength(BCHn);
	
      }
      else{
	correctedBits = noisy_codeword;
	// depad correctedBits
	correctedBits.SetLength(BCHn);
#ifdef VERBOSE
	std::cout << "Exceed BCH codes' Error Correction Capability!" << std::endl;
#endif
      }

      // write the error log into mcells
      save_error_log_in_mcells(mcells, errorLog, trees, H2, G2, metadata, encodedMetadata, numpadbits);

      // save the ber after EC
      berafter = CellUtil::compute_bit_error_rate(cells, BCHn);
#ifdef VERBOSE
      cout << "BER after EC: " << berafter << endl;
#endif

      // update the global bit error rate vector
      bit_error_rates[nosegscrubbed * 2] += berbefore;
      bit_error_rates[nosegscrubbed * 2 + 1] += berafter;

      nosegscrubbed++;
      nofaultybits = weight(codeword - correctedBits);
      if (nofaultybits == 0) {
	correctblkno++;
      }
      else {
	totalErrorBlocks++;
      }

      // increase time and move to the next segment (by increasing
      // the index i)
      // if this is the last segment, go back to scrub the first one   
      if (i == s) 
	i = 1;
      else
	i++;
      t+= gap;
      // drift all the cells containing payload data
      drift_all_cells(cells, gap);
      // drift all the cells containing metadata
      drift_all_cells(mcells, gap);      
    }

    // update performance counter
    num_seg_scrubbed += nosegscrubbed;
    average_cleanup_bit_num += bitsRefreshedBeforeEC; 
    average_mem_bit_num += bitsMemoized;
    update_common_global_performance_counters_bch(nofaultybits, correctblkno, totalErrorBlocks);
    free_PCM_cells(cells);
    free_PCM_cells(mcells);
    free_trees(trees);
  }













} // Algorithms
